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Abstract 

We consider cosmological dynamics in generalized modified gravity theory with 
the ROR term added to the action of the form R + R . Influence of RDR term 
to the known solutions of modified gravity is described. We show that in particular 
case of N = 3 these two non-Einstein terms are equally important on power-law 
solutions. These solutions and their stability have been studied using dynamical 
system approach. Some results for the case of N / 3 (including stability of de Sitter 
solution in the theory under investigation) have been found using other methods. 



o 

1 Introduction 



Theories of modified gravity become very popular in the beginning of our century 
and the main reason is attempts to understand accelerated expansion of the Universe 
which is proved observationally [1] |2]. Modifying left-hand side of Einstein equations 
of gravity it is in principal possible to obtain a late-time acceleration without invoking 
any new type of matter. It is also known that such theories can explain inflation (i.e. 
accelerated expansion at early stage of Universe evolution) using only gravitational sector 
of the theory [21 H], though inflation driven by the scalar field without any corrections to 
-^ . Einstein gravity is much popular now (nevertheless, inflation from modified gravity is not 

f^| \ currently ruled out by observations [5]). 

Standards way to modify Einstein gravity is to replace the curvature scalar R in the 
gravity Lagrangian by some scalar function f(R). In principal, it is possible to consider 
Lagrangians which depends on other curvature invariants, like the Kretchmann invariant 
RijkiR^ kl or Ricci square RikR lk etc. However, cosmological evolution in such theories 
often lead to anisotropic expansion [6l [7J [8] and, therefore, can not describe late stages of 
our Universe. This is the main reason why nowadays the class of f(R) theories is the most 
popular class of modified gravity, and we restrict ourselves by function of scalar curvature 
invariant in the present paper (it should be noted, however, that theories with scalar 
function f(G) of Gauss-Bonnet invariant G = R^mR^ —AR^R 1 ^ + R 2 are promising due 
to special properties of the invariant G [91 [TO]). 

Apart from considering scalar function of R it is possible to introduce terms depending 
on derivatives of the scalar curvature. These theories are known from the end 80th, in 



particular, theories with terms which are polynomial in the combination RDR have been 
intensively investigated [TTj [T2T IT5| 1131 [TBI fTB] . 

The goal of the present paper is to consider a theory which contains both f(R) and 
RDR modifications of the Einstein gravity. We concentrated our attention on the liner 
with respect to RDR term case, because as it was shown only liner case is the ghost free 
one [T7] . We apply theory of dynamical system to study asymptotic regimes and their 
stability. This method is used actively for finding asymptotic solutions in f(R) gravity 
and studying their stability [HI [191 ED], we extend it for theories with RDR term. We 
show that this approach works well in the situation when corrections from these two ways 
of modifying gravity are equally important. One such case is identified (f(R) = R + (3R 3 ) 
and investigated systematically. After, we discuss some problems of the used method in 
a general situation. Finally, some modifications of known cosmo logical behavior in f(R) 
theories due to presence of RDR term are described for a broader class of functions f{R). 

2 Main equations 

We consider Lagrangian in the form 

L = V=gf{R, RaR = A) + L M . (1) 

Equations of motion for such kind of theories have the form 

-\fgik + InRik - VjVfc/ij + g ik Df R + 

+f A DRR lk + n(f A R)R ik - V t V k (f A DR + D(f A R)) + D(f A DR + D(f A R))g ik + (2) 

+\Vi{f A R)V l Rg ik - Vi{f A R)VjR + \f A RURg lk = 8nGT lk , 

and its trace 

-2/ + f R R + 3Df R + Rf A DR + RD(f A R)+ 

(3) 
+3D(f A nR + D(f A R)) + V t (f A R)V l R + 2f A RDR = SttGT, 

where fn = j^, f A = qj, V«-R = J-p-. Since we interested in the flat Friedmann- 
Robertson- Walker metrics 

ds 2 = dt 2 -a 2 (t)dl 2 } (4) 

the following simple relations can be written down 

T = p - 3p, Too = P, R= -6(2H 2 + H), R 00 = -R + 3H 2 , DR = 3HR + R. (5) 
Denoting terms depending on R and RDR as f\ and fi respectively, we get 

f(R, RDR = A) = f^R) + f 2 (RDR) (6) 



As it was already noted in the present paper we study the liner case of f 2 (ROR) = aRDR, 
the function f\(R) will be specified later. It is useful to introduce the following notations: 
dfi/dR = F and a\3R = $. As F and $ are scalars, in analog with eq. fl5]) we get: 



-/i + -RF + 3H 2 F + 3HF + 2{-R + 3H 2 )<5> + 6H$ - -R z 

£ £ £ £ 



DF = 3HF + F, □$ = 3H$ + $. 
Taking into account this relations (00) component of eq. ([2]) may by written as 

Kp, 

z, z, 

and its trace ([3]) take the form 

- 2/i + FR + 9HF + 3F + 2$i? + 18#$ + 6$ + aR 2 = Kp(l - 3w), 

where we denote K = 8nG and p = wp. 

For further investigations we introduce the following dimensionless variables 



(7) 



Jb — 


6FH 2 


y = 


R 
6H 2 ' 


z = 


F 
FH 


m = 


2$ 

F ' 


n = 


2$ 
FH' 



aR 2 



n 



6FH 2 > 

Kp 
3FH 2 ' 



Dividing eq. OH]) by 3FH 2 , we rewrite it in our new variables ( ITUj) : 

x + y + z+l + m(y + 1) + n + r = Q. 
Similarly, eq. ($§ can be rewritten as 

Ax + 1y + 3z + -^ + -^§2 + 2my + 3n - 2r = £1(1 - 3w). 



(9) 



(10) 



'ir 



(12) 



Taking derivative of the introduced variables with respect to time we get the system 
of equations 



A — -R _ -rr — 9t— 

i? 6iT 3 x ^ -^#2, 



± — A o„ A 



F _ r 2 _ r H 



# FA" 2 * ^_ff 2 ' 

| = ra - mZ; (13) 

A — _?j „ 7 _ „_£. 

H ~ FH 2 nZ ,l H 2 ' 
r aRR _ _ o_ i? 



__ rtr> iyt Jj 1 

H ~ -iFH'i ' A ^ H 2 ' 



£L — Ji£ Ci 7 _ 20 — 



To complete transformation we need to express all derivatives in ([TBI such as H, R and 
other through new variables ffTUl) . First of all by using the continuity equation 

p + 3H (p + p) = (14) 

we find 3F ^ 3 = —30(1 + w). Moreover since R = —6(2H 2 + H), it easy to find that 
jp = —2 — y. Let us now introduce dimensionless parameters, characterizing the function 
h(R): 

°-^--FR' W 

„ __ RflRRR 
flRR 

Using these parameters we easily find 

R R F FF yz 

QIP ~ 6IPRF ; ^ r ~ T' 

Using (0) and ([7j) it is possible to calculate 

F z 2 c z 2 ym 



FH 2 b 2br 

Using expression (fl2|) we find 



3z. 



2$ z 2 c z 2 ym _, 

H J 4x - 2y - 2ym - 3n + 2r + 0(1 - 3w). 



FH 2 b 2br 

And finally, using R = — — 3HR we find 

aRR yzm 

~3FH 3 = V 

4 



6r. 



Note now that using eq. ( ITTj) it is possible to exclude Q from resulting system. And 
summing all above expressions and going to the conformal time, we find resulting system: 

x' = -f + x(2y + 4-z), 

y'=yi + 2y(2 + y), 

z> = ^-^-z* + z(y-l) } 

wl = n — mz, (16) 

ri = -& + ^-4x-2y-2ym+ 

+2r + (x + y + z + l+ m{y + 1) +n + r)(l — 3w) + n(y — z — 1), 

r' = -^ + r(2y-2-z), 
where we denote ' = , n d > . For some functions A there are additional relations between 

d(ln a) J x 

x and y, a well as between z and r. In the next section we consider an example for which 
such relations exist. 

3 The case of fi = R + /3i? 3 . 

In this special case some simplification can be done. First of all we note that the 
parameter c = 1. Moreover, the parameter b is a combination of variables x and y: 

b =3 x -±y. (i7) 

y 

It is possible also to find relation between r and z in this special case: 

o.z 2 y . . 



108/3(x + y) 
So we can substitute this relation into f JT6l) and exclude r from the system: 



y z 

3(x+y) 



+ 2y(2 + y), 



3(x+y) 



v +18 sm- z ^ + z ( y -ii 



(19) 



m = n — mz, 



n' = -sg^-l8^-Ax-2y-2ym-^^, 



2. 



+ (x + y + z + l + m(y + l) + n- I5 g^y )(l - 3w) + n(y - z - 1). 
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As we have two coupling constants, a and (3 our strategy is to find stationary points and 
corresponding eigenvalues separately for a = 1 and a = — 1 assuming (3 = 1. Having this 
information we then consider the general case. 

We start with a = 1. Solving the system ( 1T9|) with vanishing left-hand sides, we find 
the following stationary points. 
1. (x, y, z,m, n, ft): 

(1,-2,0,0,0,0) (20) 

Corresponding eigenvalues (which have been found numerically) are 

(-3 - 3w, -7.9040, 4.9040, -3.3677, 0.3677) (21) 

This is de Sitter point, and we can see that it is a saddle. 
2. 

Ill 3 (w 2 -l) 2(w 2 -l)(l+w) (5w 3 + 22w 2 -Mw-63)\ 

___ w , _ w __ -2-2w, 3{w _ 3 y 3(3 _ w) > 9(3 - w) ) ( } 



At least one of the eigenvalues is not negative 

2(1 + w) ^0. 
Other eigenvalues are 



(23) 



5/4w - l/4± 

l/4(-7w 2 - 50w + 313 + 8(-14w 4 - 2w 3 + 313w 2 - 624w + 387) 



1/2x1/2 



5/ Aw - l/4± 

1 /9 

l/4(-7w 2 - 50w + 313 - 8(-14w 4 - 2w 3 + 313w 2 - 62Aw + 387) 1/2 ) . 



(24) 



This is a matter dominated point in a high-curvature regime. This point is also a saddle. 
Scale factor a has the following dependence on time: 



a(t) = ao\t — to| u,+1 . 



(25) 



3. 



(0.1623, -0.4870, -6.0520, -1.0625, 6.4305, 0) (26) 

Eigenvalues are (the last eigenvalue depends on equation of state parameter w, so we 
present values found for several particular w) 

( 9.0780, w = -1 
7.5780, w = -0.5 

6.0780, w = ) (27) 

4.5780, w = 0.5 
3.0780, w = 1 



(6.4850,-1.9200,4.5650,6.0520, < 



We can see that this point is a saddle. The scale factor is 

a(t) = a \t - t | i^siao . 



4. 



;i.3732, -4.1195, 8.4778, -1.0700, -9.0710, 0) 



(28) 
(29) 



Eigenvalues are 



-13.5973, -8.0759, -8.4778, -5.5213, < 



( -12.7167, w= -1 
14.2167, w= -0.5 
15.7167,w; = 
17.2167,w = 0.5 
18.7167, w= 1 



(30) 



We can see that for all studied values of w this point is a stable node. The scale factor 
shows a phantom behavior. 

a(t) = a \t- t 1 ^T95. (31) 

5. 

(0.6979, -2.0936, 0.3742, -0.0326, -0.0122, 0) (32) 

Eigenvalues are 



(4.7703, -0.3742, -8.2381, -0.3742, < 



-0.5613, w= -1 
-2.0613, w= -0.5 
-3.4678, u> = 
-5.0613, w = 0.5 
-6.5613, w = 1 



(33) 



This is a phantom saddle with 



a(t) = a \t — t \ -"■' 



(34) 



We now turn to the case of a = — 1. Similarly as before, we get the following stationary 

points: 

1. 

(1,-2,0,0,0,0) (35) 

with eigenvalues 



-3 - 3w, -3.4198, 0.4198, -1.5000 - 5.5844i, -1.5000 + 5.5844i) 



(36) 



This is de Sitter point which is a saddle. 
2. 

1113 l-w 2 2(w 2 -l)(l+w) (5w 3 -8w 2 + 20w+A5)\ 

W,—W ,— 2 — 2w,—. r, : : , : : (37) 

2 6 '2 2' '3(w-3)' 3(w-3) ' 9{w - 3) ' 



At least one of eigenvalues is not negative one 

2(1 + w) ^0, 



(3* 
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others are 

5/ Aw - l/4± 

l/4(-7w 2 + Mw - 119 + 8(-14u? 4 + 34u? 3 + 25w 2 - 300w + 1035) 

5/4u? - l/4± 

l/4(-7u> 2 + 94w - 119 - 8(-14w 4 + 34u? 3 + 25w 2 - 300w + 1035) 

This is a matter dominated point which is a saddle with 

2 

a(t) = a \t — t Q \ w+1 . 
3. 



1/2 1/2 



1/2,1/2 



(0.7023, -2.1069, 0.4275, 0.0374, 0.0160, 0) 



(39) 



(40) 
(41) 



Eigenvalues are 



•1.7672 + 5.7921i, -0.4275, -3.5343, -1.7672 - 5.792H, { 



-0.6412, w = -1 
-2.1412, w = -0.5 
-3.5343,w = 
-5.1412, w = 0.5 
-6.6412, w= 1 



(42) 



For all studied values of w this point is an attractive focus. The behavior of the scale 
factor is of phantom type: 

a(t) = a \t - tol^^ . (43) 

All these points are summarize in the following two tables. 



Table 1. Stationary points for j\ = R + f3R 3 with /3 — 1, /2 = aR\3R with a — 1. 



Number 
of point 


Coordinates of stationary 
points 


Stability 
type 


Time dependence of 
scale factor a(t) 


1. 


x — 1, 

2/ = -2, 

2 = 0, 

m = 0, 
n = 0, 
ft = 0. 


Saddle 


a(t) = a e^i2" 


2. 


x = 1.3732, 
y = -4.1195, 

z = 8.4778, 
m = -1.0700, 
n = -9.07100, 

n = o. 


Attractive node 


U\ U 4- 1-0.4718 

a(t) = a \t — t \ 


3. 


x = 0.1623, 
y = -0.4870, 
z = -6.0520, 
m = -1.0625, 

n = 6.4305, 
tt = 0. 


Saddle 


/,\ i. , i0.6609 

a[t) = a,o\t — to 


4. 


x = 0.6979, 
y = -2.0936, 

z = 0.3742, 

m = -0.0326, 

n= -0.0122, 

tt = 0. 


Saddle 


/,\ i, , 1-10.6724 

a(t) = Oo|t — to 


5. 


X = 7j — gU>, 

y = -§ + |«;, 
z = -2 - 2iu, 

3 — 3+w ' 
_ 2 (-l+«> 2 )(l+«>) 
*" 3 -3+io ' 
n _. 1 5to 3 +22«j 2 -34w-63 
9 uj— 3 


Saddle 


2 

a(t) = ao\t — to\ w+1 



Table 2. Stationary points for f\ = R + PR 3 , (3 =, f'2 = aROR with a = —1. 



Number 


Coordinates of stationary 


Stability 


Time dependence of 


of point 


points 


type 


scale factor a(t) 




x = 1, 








y = -2, 








z = 0, 


Saddle 




1. 


171 = 0, 
71 = 0, 

ft = 0. 




a(t) = a e v^ 2 




x = 0.7023, 








y = -2.1069, 








2 = 0.4275, 


Attractive node 




2. 


m = 0.0374, 

n = 0.0160, 

ft = 0. 




/,\ 1, , 1-9.3545 
a{t) = ao\t — lq\ 




x = 2 — gW, 








2/ = ~f + ^, 








2 = -2 - 2w, 


Saddle 




3. 


= 1 -!+» 2 
'" 3 -3+w ' 

_ 2 (-l+w 2 )(l+w) 

L 3 -3+w ' 
O — 1 5w A -8w 2 +20w+45 
9 «>— 3 




2 
a(t) = a \t — t \ m+1 



It is important to compare these results with those found in the absence of f'2- term. 
For the function j\ = R + R 3 these are three stationary points - de Sitter (which is an 
unstable node), a stable vacuum power-law solution, a solution with a ~ i 1 / 2 behavior and 
an unstable matter-dominated power-law solution (the latter case represents a situation 
where the power index in not changed due to the dalambertian term). 

In the presence of the dalambertian term we have analogs of three of these four solution. 
It is interesting, that substituting in the general form of equations of motion (16) the 
function f ± = R 3 and using f 1R = 3R 2 ,f 1RR = 6R, f 1RRR = 6,6 = 2,c = 1, and, 
correspondingly, y = —3x,r = —z 2 a/72, which gives the equation of motion in the form 
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x' = ^p + x(-6x + 4 - z) 

z' = -4 - 54^ - ^(1 + 3s) 

m' = n — mz (44) 



n' = -^ + 54^f + 2a; + 6a;m - ^ + 



+(-2s + z + 1 + m(-3s + 1) + n - ^)(1 - 3w) - n(l + 3a; + 2). 

we get a bit different set of stationary points. As this theory can be treated as high- 
curvature limit of the R + R 3 theory when Einstein contribution can be neglected, we do 
not find de Sitter point. However, this is not the only difference. There are two additional 
solutions which have the following form for a = 1 for vector (x, z, m, n, fi): 

(0,0,-1,0,0) (45) 



(4,1-3^,-1,-1) (46) 



19 19 
(0,-2,--,-,0) (47) 



with eigenvalues 
and 

with eigenvalues 

(3,3-3^,-1,-1) (48) 

which correspond to a ~ t 1 / 2 behavior. Both solutions have s = 0,y = 0, and, so, they 
can not be found between solutions of Eqs. (19) due to the factor x + y in denominator. 
This is the same reason why we do not find a solution corresponding to a Friedmann 
Universe in low-curvature limit of the theory R + R 3 — this is already remarked in . 

For other values of a coordinates of one of these point changes, though the behavior 
a ~ t 1 / 2 and the eigenvalues remain unchanged which indicates that these two points are 
saddles. 

Returning to other points, we can see that the difference between a = 1 and a = — 1 
cases is in number of power-law vacuum solutions - there are three of them for a = 1 
(with one stable solution) and one (which is stable) for a = — 1. We remind a reader that 
for a = we have one stable node a ~ (t — £o)~ 10 . Using power-law ansatz we can study 
a general case. 

Substituting a(t) = a \t — t \^ into eq. (9), we get for f 1 = R + (3R 3 

F (/3(-6) 3 (2 - A) 3 + 3(A 4 (t - t ) 4 + 3/3(-6) 2 (2 - A) 2 ) + 



A<>(t-to) 



+18/3(-6)12(2 - A) 2 A + a(36(-6)(l - A) 2 (2 - A)A+ (49) 

r\ 3(l + lf) . n 

+6«(-144)(2 - A){\ - A)A 2 - 72a{2 - A) 2 A 2 - Kp A 6 {t - t ) ~ +& ) = 0. 
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Neglecting matter and Einstein term we get the equation for A 

20aA 3 + 2A 2 (15/3 - 8a) - ,4(57/3 + 6a) -6/3 = 0. 



(50) 



The solution is plotted in Fig.l. Stability analysis shows that one solution is stable, stable 
branches are shown in bold. 

Other solutions do not show any significant dependence on the value of a. The coordi- 
nates of de Sitter point do not depend on a, and stability analysis shows that it is always 
a saddle. For the matter dominated point one eigenvalue is equal to 2 — 1w for any a, so 
this point is unstable (though dimension of unstable manyfold for this point does depend 
on a). Finally eigenvalues of the points with a ~ t 1 ^ 2 does not depend on a. 





Figure 1: Power index A of a vacuum power-law solution depending on coupling constants 
a and j3. Stable branches are bold lines. 



4 The case of j x = R + R N , N ^ 3 

The equation (49) shows also why the case of f\ = R + R 3 is exceptional for theories 
with f 2 = ROR. In this case for power- law ansatz a ~ (t — to) 1 ^ 4 terms originating 
from fi and f 2 are equally important. For functions j\ = R + R N with iV > 3 in the 
high-curvature regime terms from f\ dominates (in low-curvature regime dominant terms 
are those from Einstein part of /i). This means that influence of f 2 could be important 
in some transient regime of intermediate energy, and, so, it can not be found by studies 
of stable points and corresponding asymptotic regimes. 

On the other hand, for 1 < iV < 3 the evolution of the Universe begins with R\3R 
dominant regime, after that R N corrections can be important in an intermediate regime, 
and for low-curvature stage Einstein term becomes dominant. To study such type of a 
transient dynamics we need other methods. 
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In what follows we consider two particular problems which can be solved using a 
specially chosen anzatz. 

4.1 Stability of de Sitter solution 

It is known that in a modified gravity with /i = R + R N de Sitter solution exists and 
stable for < N < 2, N =^ 1. We study what happens if we include the dalambertian 
term. 

In this section it will be better for us to use the different signature from previous one. 
Thus we put 

g ik = diag(-l, a 2 , a 2 , a 2 ) (51) 

In this notation some values change the sign, for instance we have R = Q(H + 2H 2 ), 
R 00 = 3(H + H 2 ) and respectively UR=-R- 3HR. 
We take function / in the form: 

f = R + [3R N + aRUR. (52) 

Thus 00 component of (J2J) takes the form: 

QH 2 + /3[(1 - N)R N + QH 2 NR N - 1 + QHN(N - 1)R N - 2 R] 

(53) 
+a[2RR + 3QH 3 R - R 2 - 48H 2 R - 12HR] = 0. 

This representation is more comfortable for following investigations. There is a dS-solution 
in the equation (1531 : 

<" = m^y < 54 > 

Note that here we have Rq = 12Hq. We can see that there is no dS-solution for N = 1 in 
empty space, as well as there is no dS-solution for iV = 2 in empty space. For iV = 3 we 
have dS-solution Hq = 1/(12-*//?)- Also its need to note that dS-solution exist for positive 
/3 when N > 2 and for negative one when N < 2. This result is the same as for usual 
/(i?) -gravity, because additional term aRDR can not produce new dS-solution since it 
contain only curvature derivative, however it may change the stability of dS-solutions due 
to terms with higher time derivatives. 

For finding stability conditions we study small perturbations near dS-solution Hq. 
We note that H = H + 6H, R = R + 5R = R + Q(5H + AH 5H) and R N = 
Rq + 6NR ~ 1 (5H + AHq5H), where we keeping only linear with respect to 5H terms. 
Substituting this relations into (153]) and keeping only linear terms with respect to 5H 
terms and using (|54p we find 

24H aSH^ + l2aR 5H + (lOH R a - ^^Hq 1 ) 5H 

(55) 
- (2R 2 a + 3^^) SH + 4H {N - l)8H = 0. 

We search for solutions of this equation taken in the form 8H = e xt . After substituting 
this representation into equation (155]) it transforms to usual algebraical equation of forth 
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order with four solution Aj (in the general complex case). It is clear that for stability 
solution H need -Re(Aj) < for any i. This algebraical equation is rather complicate 
one, so we use a special method (Routh-Hurwitz theorem), which detects if any root of 
this equation have a positive real part. For details see, for instance [22]. Briefly it works 
as follows: from coefficients of equation ( 155]) we construct some determinants Tj and if 
all them are positive then all real part of solutions Aj are negative. Here we present the 
results of calculations of Tf. 

T = 24# a, 



Ti 



lUH^a, 



72Hn 



4„2 



12 • 28H*a 



a- 



N{N 



N-2 



T 3 = 9 ■ 2AH I -12 2 ■ 14 • 16a 3 ^ 



**Bi^ 2 



(UN- 16) + a 



N 2 (N-iy 

(N-2) 2 



T 4 = 4H (N - 1)T 3 



(56) 

(57) 

(58) 

(59) 
(60) 



Note also that for positive a we need condition T > 0, but for negative one we need 
T < 0. First of all we can see that for N < 1 T 3 and T 4 always have opposite sing and 
therefor dS-solution is unstable for any N < 1 (note that we interested in only H q > 
case). However, for sufficiently small values of Hq all terms with high order Hq can be 
neglected and all T have the same sign in the range 1 < N < 2. It mean that in this 
range dS-solution may be stable for sufficiently big |/3|. This analytical result is confirmed 
by numerical investigation (see Fig.2). 



-N-1.2 
-N=1.6 



10 10 10 10 10 



Figure 2: Zones of stability of de Sitter solution for two values of N. The de Sitter solution 
is stable below the curves. 
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4.2 Ruzmaikin solution and its analogs 

The method we used have another limitation. Usually, having the variables determined 
we obtain the time dependence of the scale factor via Eqs. (10). The simplest way is 
to use equation for y. If the scale factor behavior is of power- law type, this equation 
gives us the exact solution. However, for other types of dynamics we can not be sure 
that the set of variables can distinguish asymptotic regimes, which are different from the 
cosmological point of view. An interesting example is the situation in modified gravity 
with /i = R + R 2 . This particular case is an exceptional one in the sense that there are no 
de Sitter solution in this theory (see above). Instead, the Ruzmaikin solution H ~ t exists 
[2~T] . From the cosmological perspective, this solution describing super-inflating Universe 
is very different from the de Sitter solution, however, it gives the same values of x, y, z as 
the latter one. Only by using initial equations (8)-(9) we can understand that something 
special happens for this particular function f\. This explains why this solution have not 
been distinguished in papers j20l [18]. It is interesting that another set of variables used 
in [6] does discriminate between de Sitter and Ruzmaikin solutions. 

In what follows we show that similar solution exists for the theory with fi = R + R 2 
and / 2 = RDR. Consider 

Bt n+1 

a(t) = const ■ e n + 1 , /gi \ 

H = Bt n . 

we get 

H = Bnt n -\ (62) 

R = -6(2B 2 t 2n + Brit* 1 - 1 ), (63) 

R = -6(4nS 2 t 2n - 1 + Bn{n - l)t n ~ 2 ), (64) 

R = -Q(4n(2n - l)SV (n_1) + Bn(n - l)(n - 2)£ n ~ 3 ), (65) 

OR = -72nBH 3n ~ 2 - 6n(lln - 7)B 2 t 2( - n ~ 1) - 6Bn(n - l)(n - 2)t n ~ 3 , (66) 

^^ = -72n(3n-2)Bh 3 ^ 1) -12n(lln-7)(n-l)B 2 t 2n ' 3 -6Bn(n-l)(n-2)(n-3)t n - 4 . (67) 

Equations of motion in the absence of matter take the form 

Et 311 - 1 + Mt 2{ - n ~^ + St 2n + 

(68) 
+D 1 t 5n - 1 + E 1 t 4n ~ 2 + S 1 t^ n ~^ + Mxt 2 ^- 2 ) = 0, 

where constant coefficients depend on coupling constants. We need two of them 

S = 3B 2 , 

(69) 
£>i = 36 ■ 12aB 5 n. 
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We denote the terms originating from the dalambertian by variables with subscripts. The 

biggest power indexes (which we should equate) belong to the terms St 2n and .D^ 5 ™ -1 . 

This means that 

2n — 5n — 1, 

3n = 1, (70) 

n = \, 

and we get the solution of Ruzmaikin type: 



4 
3B±3 



a(t) = const • e * , rj^\ 

H = BU. 

As the original Ruzmaikin solution, this solution corresponds to (x,y,z,m,n,r): 

(1,-2,0,0,0,0) (72) 

and is undistinguishable from de Sitter point from the viewpoint of our set of dimensionless 
variables. 

As the coefficient D\ comes from the dalambertian contribution, this solution exists 
also in a theory with fi — R, f 2 — RCR, i.e. theory with only dalambertian correction 
to Einstein gravity. 
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As for power-law solutions, substituting a(t) = ao(t — 1 ) 3 we can obtain solutions 
valid in the regime when R\3R dominates. There are vacuum solutions with A\ = —0.2782 
and A 2 =1.0782 (these solutions represent solutions of Eq.(50) in the limit a — > 0), and 
matter-dominated solution A = |(1 + w). These three solutions are valid in the high- 
energy regime for f\ = R + R N , 1 < N < 3. 

5 Conclusion 

We have considered cosmological dynamics in the modified gravity with both f(R) 
and R\3R terms in the action. In this theory the case of f(R) = R + R 3 appears to be 
a special one. In this particular case terms in the equations of motion originating from 
these two different modifications of Einstein gravity (namely, adding R 3 and R\3R to the 
action) are equally important on power-law solutions. This allows us to obtain several 
exact solutions and study their stability using dynamical systems approach. Our results 
show that the general picture does not change much in comparison with the known case 
of R + R 3 gravity - only one stationary point (which corresponds to a phantom solution) 
is stable, others (their number can vary) are unstable. 

Other power-law functions f(R) = R + R N lead to zones of f(R) dominance (near a 
cosmological singularity for iV > 3), or ROR dominance (R N corrections are sub-dominant 
in comparison with those from the dalambertian near a cosmological singularity if N < 3), 
as well as existence of transient epoch when some terms can be temporarily important 
(R N corrections for N < 3 and R\3R corrections for N > 3). 

However, this does not mean that sub-dominant terms in some particular epoch are 
absolutely unimportant. The reason of this is extra derivatives in terms originating from 
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RE3R contribution. This means that some stable points existing in R + R N gravity can 
loose their stability due to RDR contribution. We study in detail one particular example 
of this situation, namely stability of de Sitter solution. It is known that de Sitter solution 
exists and stable in f(R) gravity with f — R+ R N , N < 2 (JV 7^ 1). ROR terms vanishes 
on de Sitter, however, extra dimensions of the phase space may violate stability. We have 
shown that, indeed, for N < 1 the de Sitter solution is unstable (for 1 < N < 2 the 
stability depends on coupling constants). 

Finally, a modification of Ruzmaikin solution for R + R 2 gravity valid for theories with 
RE3R corrections is found. 
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